/****************************************************************************************
*                                                                                      *
*                          Projeto BioInfo - Lincs and CIn / UFPE                      *
*                                    07/10/2014                                        *
*																					   *
****************************************************************************************
* Descricao: Arquivo que contem kernel da classe Ocl.                    *
****************************************************************************************
* Responsaveis: Jefferson Ramos(jrla)													   *
*                                                                                      */
/*                                                                                     *
****************************************************************************************/


#define CUSTO_RE -1
#define CUSTO_IN -1
#define fi(ini,fim) for(i = ini; i < fim; i++)
#define fj(ini,fim) for(j = ini; j < fim; j++)
#define Byte_By_TwoBit 4
#define TAM_SEQUENCE 240
#define MASK_NCL1 0x03
#define MASK_NCL2 0x0C
#define MASK_NCL3 0x30
#define MASK_NCL4 0xC0
#define CUSTO_TRACEBACK_MATCH 3
#define CUSTO_TRACEBACK_MISSMATCH -2
#define CUSTO_TRACEBACK_GAP_OPEN -5
#define CUSTO_TRACEBACK_GAP_EXTENSION -2
#define PE_UNROLL()\
		score_diag2 = score_up[0];\
		score_diag3 = score_up[1];\
		score_diag4 = score_up[2];\
		score_left = score_vector[i + score_vector_ID];\
		score_vector[i - 1 + score_vector_ID] = score_up[3];\
		score_result = maximum(	score_diag1 + custo_diag[0], score_diag2 + custo_diag[1],\
								score_diag3 + custo_diag[2], score_diag4 + custo_diag[3],\
								custo_diag, score_up, score_left,\
								&store_direction);\
		score_diag1 = score_vector[i + score_vector_ID];\
		score_vector[i + score_vector_ID] = score_result;\
		input_direction[(index_byte-1) + (i-1)*(TAM_SEQUENCE/Byte_By_TwoBit)]= store_direction;

#define __FAST_RELAXED_MATH__

int maximum(int diag1, int diag2, int diag3, int diag4,
			int *cust_diag,int *score_up, int left,
			unsigned char *direction){
			
	char direction_local;
	if (diag1  < (left + CUSTO_IN)) //Mach ganha de insersao
	{
		if ((left + CUSTO_IN) >= (score_up[0] + CUSTO_RE)) {
			direction_local = 1;
			score_up[0] = (left + CUSTO_IN); //Insersao
		}
		else {
			direction_local = 0;
			score_up[0] = (score_up[0] + CUSTO_RE); //Remocao
		}
	}
	else //Insersao ganha do Mach
	{
		if (diag1 < (score_up[0] + CUSTO_RE)) {
			direction_local = 0;
			score_up[0] = (score_up[0] + CUSTO_RE); //remocao
		}
		else {
			if (cust_diag[0] == 1) {
				direction_local = 3;
				score_up[0] = diag1; //mach
			}
			else {
				direction_local = 2;
				score_up[0] = diag1; //mach
			}
		}
	}
	left = score_up[0];

	if (diag2  < (left + CUSTO_IN)) //Mach ganha de insersao
	{
		if ((left + CUSTO_IN) >= (score_up[1] + CUSTO_RE)) {
			direction_local |= 1<<2;
			score_up[1] = (left + CUSTO_IN); //Insersao
		}
		else {
			direction_local |= 0<<2;
			score_up[1] = (score_up[1] + CUSTO_RE); //Remocao
		}
	}
	else //Insersao ganha do Mach
	{
		if (diag2 < (score_up[1] + CUSTO_RE)) {
			direction_local |= 0<<2;
			score_up[1] = (score_up[1] + CUSTO_RE); //remocao
		}
		else {
			if (cust_diag[1] == 1) {
				direction_local |= 3<<2;
				score_up[1] = diag2; //mach
			}
			else {
				direction_local |= 2<<2;
				score_up[1] = diag2; //mach
			}
		}
	}
	left = score_up[1];

	if (diag3  < (left + CUSTO_IN)) //Mach ganha de insersao
	{
		if ((left + CUSTO_IN) >= (score_up[2] + CUSTO_RE)) {
			direction_local |= 1<<4;
			score_up[2] = (left + CUSTO_IN); //Insersao
		}
		else {
			direction_local |= 0<<4;
			score_up[2] = (score_up[2] + CUSTO_RE); //Remocao
		}
	}
	else //Insersao ganha do Mach
	{
		if (diag3 < (score_up[2] + CUSTO_RE)) {
			direction_local |= 0<<4;
			score_up[2] = (score_up[2] + CUSTO_RE); //remocao
		}
		else {
			if (cust_diag[2] == 1) {
				direction_local |= 3<<4;
				score_up[2] = diag3; //mach
			}
			else {
				direction_local |= 2<<4;
				score_up[2] = diag3; //mach
			}
		}
	}
	left = score_up[2];

	if (diag4  < (left + CUSTO_IN)) //Mach ganha de insersao
	{
		if ((left + CUSTO_IN) >= (score_up[3] + CUSTO_RE)) {
			direction_local |= 1<<6;
			score_up[3] = (left + CUSTO_IN); //Insersao
		}
		else {
			direction_local |= 0<<6;
			score_up[3] = (score_up[3] + CUSTO_RE); //Remocao
		}
	}
	else //Insersao ganha do Mach
	{
		if (diag4 < (score_up[3] + CUSTO_RE)) {
			direction_local |= 0<<6;
			score_up[3] = (score_up[3] + CUSTO_RE); //remocao
		}
		else {
			if (cust_diag[3] == 1) {
				direction_local |= 3<<6;
				score_up[3] = diag4; //mach
			}
			else {
				direction_local |= 2<<6;
				score_up[3] = diag4; //mach
			}
		}
	}
	*direction = direction_local;
	return score_up[3];
}

__kernel void
nw_kernel1(	__global char * query,
			__global  int * score_vector,
			__global char * bank,
			__global  int * output_taxa,
			const int alignment_by_individuals,
			const int num_individuals
)
{
	unsigned char input_direction[TAM_SEQUENCE*TAM_SEQUENCE/Byte_By_TwoBit];//TAM_SEQUENCE
	
	size_t global_ID = get_global_id(0);
	
	int query_ID = (global_ID%alignment_by_individuals) * (TAM_SEQUENCE);
	int score_vector_ID = (global_ID%(alignment_by_individuals*num_individuals)) * (TAM_SEQUENCE+1);
	
	int score_diag1, score_diag2, score_diag3, score_diag4;
	int score_left, score_result, score_cont;
	int custo_diag[4];
	int score_up[4];
	char individuals_nucleotideo[4];
	
	int j,i;
	int index_byte;
	
	int direction;
	int direction_anterior;
	int taxa;
	
	char query_nucleotideo;
	unsigned char store_direction;
	
	for(j = 0; j < TAM_SEQUENCE; j++){
		score_vector[j + score_vector_ID] = CUSTO_IN*(j % (TAM_SEQUENCE+1));
//		query_l[j] = query[query_ID + j];
	}
	score_vector[j + score_vector_ID] = CUSTO_IN*(j % (TAM_SEQUENCE+1));

	for(index_byte = 1, score_cont = 0, j = 1; j <= TAM_SEQUENCE; j += 4, score_cont -= 4, index_byte++)
	{
		score_up[0] = score_cont + CUSTO_RE;
		score_up[1] = score_up[0] + CUSTO_RE;
		score_up[2] = score_up[1] + CUSTO_RE;
		score_up[3] = score_up[2] + CUSTO_RE;

		individuals_nucleotideo[0] = bank[(global_ID*TAM_SEQUENCE) + j - 1];
		individuals_nucleotideo[1] = bank[(global_ID*TAM_SEQUENCE) + j];
		individuals_nucleotideo[2] = bank[(global_ID*TAM_SEQUENCE) + j + 1];
		individuals_nucleotideo[3] = bank[(global_ID*TAM_SEQUENCE) + j + 2];

		score_diag1 = score_vector[score_vector_ID];
		
		for(i = 1; i <= TAM_SEQUENCE; i++){
			query_nucleotideo = query[i-1 + query_ID];

			custo_diag[0] = (individuals_nucleotideo[0] ^ query_nucleotideo) ? -1 : 1;
			custo_diag[1] = (individuals_nucleotideo[1] ^ query_nucleotideo) ? -1 : 1;
			custo_diag[2] = (individuals_nucleotideo[2] ^ query_nucleotideo) ? -1 : 1;
			custo_diag[3] = (individuals_nucleotideo[3] ^ query_nucleotideo) ? -1 : 1;
			PE_UNROLL();
		}
	}

	i = TAM_SEQUENCE-1;
	j = (TAM_SEQUENCE/Byte_By_TwoBit)-1;
	direction_anterior = 3;
	direction = 3;
	taxa = 0;

	short status = 5;
	bool valid = 1;
	
	while (valid) {	//traceback
		if(direction != 1 || status == 1)
			store_direction = input_direction[(j) + (i)*(TAM_SEQUENCE/Byte_By_TwoBit) ];

		switch(status){//states
			case 4:
				status--;
				direction = ((MASK_NCL3 & store_direction)>>4);
			break;
			case 3:
				status--;
				direction = ((MASK_NCL2 & store_direction)>>2);
			break;
			case 2:
				status--;
				direction = ((MASK_NCL1 & store_direction));
			break;
			default:
				status = 4;
				direction = ((MASK_NCL4 & store_direction)>>6);
			break;
		}


		switch(direction){
			case 0://ups
				if (i - 1 >= 0){
					i = i - 1;
					status++;
					if(direction_anterior == 1 || direction_anterior == 0)
						taxa += CUSTO_TRACEBACK_GAP_EXTENSION;
					else
						taxa += CUSTO_TRACEBACK_GAP_OPEN;
				}
				else
					valid = 0;
			break;
			case 1://lefts
				if(status != 1){
					if(direction_anterior == 1 || direction_anterior == 0)
						taxa += CUSTO_TRACEBACK_GAP_EXTENSION;
					else
						taxa += CUSTO_TRACEBACK_GAP_OPEN;
				}
				else{
					if(j - 1 >= 0)
					{
						j = j-1;
						if(direction_anterior == 1 || direction_anterior == 0)
							taxa += CUSTO_TRACEBACK_GAP_EXTENSION;
						else
							taxa += CUSTO_TRACEBACK_GAP_OPEN;
					}
					else
						valid = 0;
				}
			break;
			default://matchs e missmatchs
				if(direction == 2)
					taxa += CUSTO_TRACEBACK_MISSMATCH;
				else
					taxa += CUSTO_TRACEBACK_MATCH;

				if(i - 1 >= 0)
				{
					if(status != 1)
						i = i - 1;
					else
					{
						if (j - 1 >= 0){
							i = i - 1;
							j = j - 1;
							status = 5;
						}
						else
							valid = 0;
					}
				}
				else
					valid = 0;
			break;
		}
		direction_anterior = direction;
	}

	status--;
	if(j > 0)
	{
		if(direction_anterior == 1 || direction_anterior == 0)
			taxa += (CUSTO_TRACEBACK_GAP_EXTENSION * j * 4) + CUSTO_TRACEBACK_GAP_EXTENSION * (status);
		else
			taxa += CUSTO_TRACEBACK_GAP_OPEN + (CUSTO_TRACEBACK_GAP_EXTENSION * j * 4) + CUSTO_TRACEBACK_GAP_EXTENSION * (status-1);
	}else if(j == 0 && status != 0){
		if(direction_anterior == 1 || direction_anterior == 0)
			taxa += (CUSTO_TRACEBACK_GAP_EXTENSION * (status));
		else
			taxa += CUSTO_TRACEBACK_GAP_OPEN + (CUSTO_TRACEBACK_GAP_EXTENSION * (status-1));
	}
	if(i > 0){//completando traceback pelas linhas (traceback termina pelo lado)
		if(direction_anterior == 1 || direction_anterior == 0)
			taxa += CUSTO_TRACEBACK_GAP_EXTENSION * i;
		else
			taxa += CUSTO_TRACEBACK_GAP_OPEN + (i-1)*CUSTO_TRACEBACK_GAP_EXTENSION;
	}
	
	output_taxa[global_ID] = taxa;

	return;
}
